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Abstract. A Bose-Hubbard Hamiltonian, modelling cold bosons in an optical 
lattice, is used to simulate the dynamics of interacting open quantum systems as 
subsystems a larger closed system, avoiding complications like the introduction 
of baths, complex absorbing potentials or absorbing boundaries. The numerically 
exact unitary dynamics is compared with effective descriptions of the subsystems 
based on non-Hermitian Hamiltonians or Lindblad master equations. The validity 
of popular models with constant decay rates is explicitly analyzed for decaying 
single and double wells. In addition we present a discrete lattice version of the 
Siegert approximation method for calculating decay rates. 
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1. Introduction 

The behaviour of interacting quantum particles in open systems is of fundamental 
interest. Such systems have been studied experimentally and theoretically in different 
contexts, including electronic transport in semiconductors and nanostructures [T] and 
cavity QED [g. 

Cold atom experiments [3J|31[51|B] permit a new way to study interacting quantum 
particles in open systems and thus a way to test different theoretical approaches 
for describing these systems, like effective non-Hermitian Hamiltonians or master 
equations. Benefits of the cold atom approach include the possibility of creating 
various trap geometries, tunability of the interaction between the particles and the 
absence of defects and impurities. In particular one can realize simple setups where 
open systems can be studied as smaller subsystems of a larger, closed system the 
dynamics of which is governed by familiar Hermitian Hamiltonians leading to unitary 
time-evolution. In the present article this concept is implemented theoretically for 
the particular situation of tunnelling decay of cold bosons within an optical lattice, 
modeled by the Bose-Hubbard Hamiltonian. 



In recent years the decay dynamics of trapped cold bosons was considered in 
various contexts. The nonexponential tunnelling decay of Bose-Einstein condensates 
was analyzed in the Gross-Pitaevskii mean-field approximation for different trap 
geometries [71 [8j HI [TOl HH [12] • Few boson tunnelling was considered in [13] by means 
of a numerically exact method yielding deviations from the mean-field behaviour 
for small particle numbers. Furthermore, the decay of bosons in optical lattices 
has been studied using effective non-Hermitian Hamiltonians [HI [TS] [M] , Lindblad 
master equations and the Bogoliubov Backreaction approximation [171 ITS] . In these 
works penomenological models with localized constant decay rates are used. While 
the latter appear appropriate for describing decay mechanisms like, e.g. particle loss 
due to a focused laser beam or electron beam (the latter was implemented successfully 
in a recent experiment [H]), this does not have to be the case for tunnelling decay 
considering the nonexponential behaviour found in some of the studies mentioned 
above. 

Here we analyze the tunnelling decay of cold bosons in a lattice within the 
framework of the onc-dimensional Bose-Hubbard model 




Figure 1. Bosc-Hubbard lattice where two sites (filled circles) forming a double 
well with internal tunnelling coefficient J are weakly coupled (coefficient oi ^ f!) 
to a long chain of sites (empty circles) with internal tunnelling coefficient Q. 
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where Uj is the annihilation operator of a particle in the lattice site j and ej are on-site 
single particle energies. The interaction parameters Uj depend on the shapes of the 
local ground states in the respective lattice sites and the local tunnelling rates Jj on 
the overlap of the ground states in the respective adjacent sites pU] . 

For a high total number of bosons N, the operators Oj can be replaced by complex 
numbers VTV Cj representing their respective coherent state expectation values (see 
e.g. [H]) which leads to the mean-field Hamiltonian 

H = J2 Un\c,\' ~ (c;c,+i + c*+ic,) + ^N{N l)\cA .(2) 
j 

The dynamics of the on-site amplitudes Cj is then given by icj = dH/dc* which yields 
the system of coupled discrete nonlinear Schrodinger equations or Gross-Pitaevskii 
equations 

iCj = e,c, - ^c,+i - %icj_i + Uj{N - l)|cjf c, (3) 

in scaled units with h = 1 that we will use throughout this article. In the following 
the dynamics obtained from numerically exact solutions of equations ^ and is 
compared with the results of effective theoretical models, including non-hcrmitean 
Hamiltonians and master equations, describing small subsystems that are weakly 
coupled to the rest of the lattice. This is illustrated in figure [1] for a double well 
subsystem. 

The main objectives of this article can be summed up as follows: 

• An experimentally realizable system, namely cold bosons in a lattice, is used 
to study non-Hermitian quantum dynamics of interacting particles under clean, 
well-controlled conditions. 

• Effective descriptions of open systems, like e.g. non-Hermitian Hamiltonians or 
Lindblad master equations, are compared with numerically exact calculations 
within genuinely closed systems without any additional approximations like 
complex absorbing potentials, absorbing boundaries or the introduction of particle 
baths. 

• The validity of popular phenomenological models with constant decay rates is 
explicitly tested for a concrete mechanism, namely tunnelling decay within a 
Bose-Hubbard lattice. 

• The full Bose-Hubbard dynamics is compared with the mean-field approximation. 

• The relative technical simplicity of our approach makes sure that the results are 
neither obscured nor compromized by mathematical or numerical subtleties. 

• Methodically, the discrete lattice version of the Siegert approximation method is 
presented as a technically simple alternative to decay rate calculations based on 
Green functions or Fermi's Golden Rule. 

In section [5] we consider single well tunnelling, i.e. tunnelling out of one site 
coupled to a long chain, double well tunnelling will be analyzed in section [31 

2. Single well tunnelling 

We consider a situation where one site, let us say site 0, is coupled weakly to a long 
chain of sites, i.e. we choose the tunnelling coefficients Jj ~ il, j > and Jq = cj <C il. 
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For the sake of simplicity we further choose ~ e, Uq ^ U and ej ~ 0, C/j = for 
J > 0. such that there is no interaction within the long chain. 

In the following we will find an approximation to the decay coefficient for 
tunnelling from the single site into the chain by means of the Siegert approximation 
method [2H [TUl which we adapt for use in a discrete lattice. 

By means of the usual ansatz Cj(t) = Cj exp(— i^<) with the chemical potential /i, 
equation Q leads to the time-independent nonlinear Schrodinger equations 

MCo = -|ci+eco + C/(7V-l)|copco (4) 

MCl = --^Co - -C2 (5) 

MCj = -^(cj-i + Cj+i) , .7 > 1 • (6) 

For J > 1 we make an outgoing wave (Siegert) ansatz Cj = Aexp(ifcj) in analogy to 
outgoing plane waves in continuous space. Thus equation ([5]) leads to the dispersion 
relation /^(fc) = --ficos(fc). In order to obtain the current in the region j > 1 we 
consider 9t|cjp = c*Cj + c*Cj (cf. e.g. [24]). From ([6]) we obtain 

|Cj|2 = -i^ (c*+lCj - C*C.3 + x + C*_^C.j - C*Cj_i) = - + Jjj-l(7) 

where we have identified the currents J^jj+i ~ iil{c*^iCj — c*Cj+i)/2 and Jj,j-i = 
iil(c*_]^Cj — c*Cj_i)/2 going from site j to sites j + 1 and j — 1 respectivelly. 
Inserting Cj = j4exp(ifcj) the "outgoing" current becomes — — l^|Apsin(fc) = 

— — (p/n)'^ =: J' where we have used the dispersion relation in the last step. 
The amplitude A follows from equation ([Sj. Iserting Cj = Aexp(ifcj), j > 1 and 
the dispersion relation leads to A = cqlu/D,. Particle number conservation requires 
I Cop = J'. If we assume an exponential decay, the decay rate must be given by 
r -|co|2/|coP = J/\co\^ which leads to 



Due to w ^ the term — cJCi/2 = — (a;^/2r2)co exp(ifc) in ^ can be neglected so 
that the chemical potential is approximately given by /i « e + U{N — l)|cop. The 
occupation of the single site thus satisfies the differential equation 



Alternatively, this result can be obtained by means of Green functions (see, e.g. |25j ) 
or Fermi's Golden Rule [26]. Equation ([9]) implies that due to the dispersion relation 
in the lattice there is no decay for |e + U{N — l)|cop| > |fi| which can be confirmed 
numerically. 

Now we turn to the full Bose-Hubbard model. We attempt an effective description 
by means of a master equation for the probabilities Pm that the single site subsystem 
is occuppied by M particles, < M < N [27]. We apply the sequential tunnelling 
approximation, i.e. we neglect the possibility of simultaneous decay of two or more 
particles, concentrating on single particle processes only [57]. As in the mean-field limit 
the transition rates can be obtained using the Siegert approximation method. Instead 
of the mean-field system ((l])-((n]) we consider the Heisenberg equations ioj = [aj,H] 
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derived from the Haniiltonian ([T]). Since there is no interaction in the chain the 
current can be obtained in complete analogy with the mean-field case. The continuity 
equation still holds so that the tunnelling rate par particle ist still given by (|8]). The 
single site system goes from an M-particle state to an M — 1 particle state when a 
particle tunnels into the chain. At site the Hcisenbcrg equation reads 

W 4-2 

iflo = --ai + eoflo + C/flo flo • (10) 

The chemical potential for going from state \M) to state |M — 1) is approximately 
determined by its stationary version 

^iM{M - l\ao\M) « (e + C/(M - 1))(A/ - l\ao\M) (11) 

where we have neglected the small first term due to a; ^ in analogy to the mean- 
field case. Thus the chemical potential reads fiM ^ e + U {M — 1). The corresponding 
tunnelling rates per particle are the given by 

^''-^V ^^^^ 

for states with M particles. This leads to the rate equation model 

PN{t) = -NTNPNit) , Po{t) - TiPiit) , (13) 

PM{t) = {AI + l)rM+iPM+i{t)-MTMPMit) , 1<M< iV-l(14) 

for the relative occupations Pm of the A/-particlc states. The total occupation of the 
first site is then given by 

N 

no{t)=Y,MPM{t). (15) 
The system can be integrated analytically 

PN{t)^cM-^Nt)PNiO), Po(t) - / A(T)dr, (16) 



PM{t)= f eM-MTMit-T))iM + l)TM+iPM+idT, 1<M<N-1.{17) 
Jo 

Now we compare the approximate descriptions obtained above with numerically 
exact integrations of the full mean-field system and the full Bose- Hubbard system for a 
finite lattice. We start with all particles in the first site. Decay without backrcfiection 
can be observed for short times when there is no backreflection at the end of the 
finite chain. The simulation time Ts in a chain of length L is thus determined by 
(yt/2)Ts ~ L, Ts ~ 2L/Vt with the phase velocity £7/2. The results are shown in 
figure m In accordance with the results obtained in [T3] for single- well tunnelling in a 
different (non-lattice) setup, there is a clear difference between few-particle tunnelling 
and mean-field tunnelling. In both cases the numerically exact results (solid lines) 
are quite well described by the respective approximations discussed above (dashed 
lines). For N = 2 particles there are some deviations, which might occur due to 
pair-tunnelling, which has been neglected in the in the rate equation model. For both 
= 3 particles and the mean-field there is a much better correspondence between the 
approximations and the exact results, which might be due to the diminished relative 
importance of pair-tunnelling. For comparison, the right panel also shows the rate 
equation prediction for = 15 particles (dashed dotted line) which almost coincides 
with the mean-field results which indicates that the mean-field limit is approached 
with relatively few particles. 
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Figure 2. (Color online) Decay from a single site with parameters e = 0, lo = 0.1, 
f2 = 1, h = 1. The occupation no(t)/no{0) of the first lattice site is shown as 
a function of time. Numerically exact results in a lattice with 200 (left panels) 
and 150 sites (right panels) (solid lines) are compared with the corresponding 
approximations (dashed lines). Left panel from top to bottom: N = 2 particles 
with U (N — 1) = 0.8; mean-field with U {N — 1) = 0.8 and interaction-free system 
with U = 0. Right panel from top to bottom: N = 3 particles with U{N—1) = 0.8; 
N = 2 particles with U{N - 1) = 0.8; mean-field with U{N - 1) = 0.8; dashed- 
dotted line: Rate equation result for = 15 particles with U{N — 1) = 0.8. 



3. Double well tunnelling 

Now we consider an open double well, consisting of two sites and 1 weakly coupled 
to a long chain of sites. The quasistationary states can again be obtained by means of 
the Siegert approximation method. The nonlinear Schrodinger equations describing 
our system in the mean-field limit now read 

Atco = -^ci+eoCo + C/(iV-l)|copco (18) 
fici = -|c2 - ^co + eici + U{N - l)\ci\^ci (19) 



n 



MC2 = -"^Ci - -C3 (20) 

ficj = + Cj+i) , j >2. (21) 

The outgoing wave solution, dispersion relation and current in the chain (j > 2) 
are still given by Cj = Aexp(ifcj), /i(fc) ~ — i7cos(fc) and ~ ~ft\A\'^ ^yl — {fi/il)"^ 
respectivelly. Analogous to the single site case the amplitude A is obtained from 
equation (PU)) as A = ciw/rt. The continuity equation |cop + |cop = for the first 
two sites then yields the double well decay rate 



r = ^ = l^il' if! (22) 

' |C0P + |C1|2 |co|2 + |ciP J] V ^ ' 

where the chemical potential ^ is approximately determined by equations (jlSp and (|19p 
the coupling to the chain is neglected, i.e. a; is set to zero. Often, phenomenological 
models with constant decay rates are used to describe open systems. For our single 
site model with tunnelling decay from the previous section this is justified for /i <C 
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when r w Yf 7- Physically this means that the internal dynamics of the subsystem 
is slow compared to the transport velocity in the chain. In this limit, one can try to 
model our decaying double well by simply adding an imaginary constant — i7/2 to the 
on-site energy of the second well which leads to the model 

ico - -^ci + eoco + U{N ~ l)|copco (23) 

ici = -|co + (ei - i^)ci + U(N - \)\cx\^c^ (24) 

considered in a number of recent works (see, e. g. [T71 [THl [H])- It is interesting to 
compare the predictions of this model for the decay rates of stationary states with the 
Siegert approximation result (j22p . To this end we treat the term — i7Ci/2 as a small 
perturbation. The unperturbed chemical potential and amplitudes co and c\ are then 
determined by the stationary states of the system ((23)) and (j24[) with 7 set to zero, 
obtained in the usual way with the ansatz Cjify = exp(— i/ii)cj, 7 = 0,1. The first 
order perturbation theory correction is straightforwardly obtained as 

This purely imaginary correction corresponds to a decay rate 

r;, = -2Im(AM) = I 1, 7 = I il'l'i (26) 

|cor + |cir |cor + kir ^ 

which is equal to ([^^ in the limit /i ^ il. Therefore both descriptions are compatible 
in this limit. This result can be straightforwardly generalized to any finite number of 
sites coupled to a chain. 

The model given by ([23| . ([24)) with 7 = is a well-studied system whose 
eigenvectors and corresponding stability properties can be obtained analytically [55]. 
In the particularly simple symmetric case with eo = = ei the eigenvalues and 
eigenvectors are given by 

/x = -Jry/2 4-C/(A^-l)™?V(l + '7^) (27) 

with 

77e{±l}, |C/(7V-l)n| < |J| (28) 

and 



r/€ {±1,-C/(A^- l)n/J± ^V^{N ^Vfn^jJ'^ - 1}, \U{N -\)n\ > |J|, (29) 

where n = |coP -I- |cip is the norm of the eigenvector and rj = cq/ci is real and 
nonzero. For a repulsive interaction U > Q the anti-symmetric solution with i] ~ —1 
becomes dynamically unstable for \U{N — \)n\ > \ J\. The two solutions that only 
exist for sufficiently interactions are strongly localized in either the left or the right 
site, thus breaking the symmetry of the system. Figure [3] shows the time propagation 
of two of the eigenstates of the symmetric double well due to the full mean-field 
dynamics of the whole system including the chain (solid lines) and according to the 
model (dashed dotted lines). The two approaches coincide well for both the 

symmetric ground state (upper panels) and the eigenstate localized in the right well 
(lower panels). For the symmetric ground we find an almost quasistationary decay 
behaviour since both the occupation |ci(t)p of the second lattice site and the total 
occupation n{t) = \co{t)\^\ -(- |ci(t)p of the double well show almost pure decay with 
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Figure 3. (Color online) Decay from a double well with parameters ei = 0e2, 
J = 0.5, uj = 0.25, n = 3, h = 1. U{N - 1) = 0.8 for the symmetric 
ground state (upper panels) and the eigenstate localized in the right well (lower 
panels). The occupation |ci(t)p of the second lattice site and the occupation 
n{t) = |co(OPl + |ci(i)P of the double well are shown as functions of time. 
Predictions of the non-Hermitian two mode model (dashed dotted lines) are 
compared with a numerically exact propagation of the full nonlinear Schrodinger 
equation in a lattice with 200 sites (solid lines). The dashed lines in the right 
panels correspond to quasistationary decay with decay rate II26II . 



only a small oscillation. Consequently the total occupation n(t) is well described by 
a decay with rate Fjj (dashed lines). For the symmetry breaking state we observe 
a quasistationary decay until the effective interaction U{N — l)n{t) drops below the 
threshold value | J| for the existence of symmetry breaking eigenstates and the site 
occupations start to oscillate. In spite of this fact the total occupation n{t) is still 
reasonably well described by a decay with rate (dashed line). This behaviour is 
in agreement with previous studies of open double well systems containing detailed 
discussions of the dynamics of the (quasi-) eigenstates [ITj [181 [13 El- Figure [4] 
shows the time propagation for the initial conditions co(0) = 0, ci(0) = 1, i.e. all 
particles are initially in the right well. For a moderate value U{N — 1) = 0.8 of the 
interaction (upper panels) the oscillatory behaviour of the full system (solid lines) 
is well capturd by the non-Hermitian two mode model ([25)) . ([^^ (dashed dotted 
lines). The lower panels display the dynamics for a higher value U{N — 1) = 1.1 
of the interaction which lies above the treshold U{N — 1) = 2\J\ = 1 for running 
phase self-trapping (see e.g. [551 [ISl 1301 [SB [321 [33) where interactions prevent a 
total population transfer between the wells in the corresponding hermitian double 
well system. Such an oscillation with a small amplitude can also be observed here for 
short times. However, after the first oscillation the effective interaction U{N — i)n{t) 
drops below the threshold value 2\J\ — 1 so that oscillations with larger amplitudes 
are possible again. Before the threshold is reached, the non-Hermitian two mode 
model (dashed dotted lines) provides an excellent description of the dynamics of the 
full system (solid lines), afterwards deviations occur but the qualitative behaviour ist 
still correctly described. 

Now we again turn to the dynamics of the full Bosc-Hubbard system and compare 
it with an effective two mode description of our open double well, namely the Lindblad 
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Figure 4. (Color online) Decay from a double well with parameters ei = 0e2, 
J = 0.5, u = 0.25, n = 3, h=l, U{N - 1) = 0.8 (upper panels), U{N - 1) = 1.1 
(lower panels) for the initial conditions co(0) = 0, ci(0) = 1. The relative 
occupation |ci(t)p of the second lattice site and the relative occupation n{t) = 
I CO (i ) I ^ I + I ci (t ) p of the double well are shown as functions of time. Predictions 
of the non-Hermitian two mode model (dashed dotted lines) are compared with 
a numerically exact propagation of the full nonlinear Schrodinger equation in a 
lattice with 200 sites (solid lines). 



master equation [33] 

P = -i[ffo, p] - ]p {a\aip + pa\ai - 2ai/5aJ^ (30) 

for the density matrix p where the second term describes constant decay from site 1 
with rate 7 whereas the first term provides the hermitian part of the time evolution 
with the two site Bose-Hubbard Hamiltonian 

Ho = ~{Jl2){a\ao + ajai) + {U/2){al^al + a\^a\) . (31) 

The time-dependent expectation values of an operator O is then given by the trace 
{0{t)) = Tr(p(t)0). Recently it was shown [T7l[T8] that the Lindblad master equation 
pop reduces to the nonlinear non-Hermitian model (P3)) . in the mean- field limit. 
For a small number of particles, equation pop can be integrated directly, for higher 
particle numbers not considered here it can be solved using Monte Carlo methods 
[53] . Figure |S] displays the relative occupation of the second lattice site Ni{t)/N, 
Ni{t) = Tr(p(t)a|ai) (left panels) and the relative total occupation {No{t) + Ni{t))/N 
of the double well (right panels) for initial conditions where only the second well 
is occupied at t = 0. Both for iV = 2 particles (upper panels) and for = 3 
particles the full Bose-Hubbard dynamics (solid lines) is well reproduced by the 
Lindblad master equation (dashed dotted lines). Compared to the mean- field case 
the oscillatory dynamics of the system appears less regular due to the occurrence 
of various frequencies corresponding to transitions between different eigenstates of 
the many- (or here rather few-) particle system, an effect which, if regarded from a 
mean-field point of view, is also referred to as quantum fluctuations in this context 

[30] HUES]- 
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Figure 5. (Color online) Decay from a single site with parameters ei = 0e2, 
LU = 0.25, Q = 3, h = 1,U{N — 1) = 0.8 for initial conditions where only the 
second well is occupied at t = 0. The relative occupation Ni(t)/N of the second 
lattice site and the relative occupation {No{t) + Ni{t))/N of the double well are 
shown as functions of time. Predictions of the Lindblad master equation I I30I I 
(dashed dotted lines) are compared with a numerically exact propagation of the 
full Bosc-Hubbard Hamiltonian in a lattice with 150 sites (solid lines). Upper 
panels: N = 2 particles, lower panels: N = 3 particles. 



4. Conclusion 

In this article, an experimentally realizable system, namely cold bosons in an optical 
lattice modeled by a Bose-Hubbard Hamiltonian was used to theoretically investigate 
the dynamics of open interacting many-particle systems within a closed system 
setup. Open single and double well systems were simulated as one respectivelly two 
sites weakly coupled to a long but finite Bose-Hubbard chain avoiding additional 
approximations due to absorbing boundaries, complex absorbing potentials or the 
introduction of baths. 

Even for single site tunnelling deviations of the mean-field dynamics from the 
full Bosc-Hubbard dynamics were found in accordance with reference |13j . The non- 
exponential decay behaviour and in particular the differences between Bose-Hubbard 
and mean-field were demonstrated to depend on the dispersion relation in the chain. 
Both in the mean-field and many-particle case the dynamics was well described by rate 
equation models derived using a discrete lattice version of the Siegert approximation 
method. 

It was shown that a description of tunnelling decay by means of constant local 
decay terms is justified if the chemical potential of the considered subsystem is small 
compared to the tunnelling coeSicient in the lattice. In the latter limit the dynamics of 
an open double well was analyzed. For the latter it was found that full mean-field and 
Bose-Hubbard dynamics is well described by a non-hermitian nonlinear Schrodinger 
equation respectivelly a Lindblad master equation with a constant decay term. 

In summary, decaying interacting quantum systems were analyzed for a concrete 
physical situation by means of a closed system approach, explicitly confirming the 
validity of popular effective theoretical descriptions that are usually applied on a 
phenomenological basis. 
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